arXiv:1503.07632v2 [math.NA] 27 Mar 2015 


WELL-CONDITIONED FRACTIONAL COLLOCATION METHODS USING 
FRACTIONAL BIRKHOFF INTERPOLATION BASIS 

YUJIAN JIAO * 1 , LI-LIAN WANG 2 * AND CAN HUANG 3 


Abstract. The purpose of this paper is twofold. Firstly, we provide explicit and compact formu¬ 
las for computing both Caputo and (modified) Riemann-Liouville (RL) fractional pseudospectral 
differentiation matrices (F-PSDMs) of any order at general Jacobi-Gauss-Lobatto (JGL) points. 
We show that in the Caputo case, it suffices to compute F-PSDM of order p E (0,1) to com¬ 
pute that of any order k + p with integer k > 0, while in the modified RL case, it is only 
necessary to evaluate a fractional integral matrix of order p E (0,1). Secondly, we introduce 
suitable fractional JGL Birkhoff interpolation problems leading to new interpolation polynomial 
basis functions with remarkable properties: (i) the matrix generated from the new basis yields 
the exact inverse of F-PSDM at “interior” JGL points; (ii) the matrix of the highest fractional 
derivative in a collocation scheme under the new basis is diagonal; and (iii) the resulted linear 
system is well-conditioned in the Caputo case, while in the modified RL case, the eigenvalues of 
the coefficient matrix are highly concentrated. In both cases, the linear systems of the collocation 
schemes using the new basis can solved by an iterative solver within a few iterations. Notably, 
the inverse can be computed in a very stable manner, so this offers optimal preconditioners for 
usual fractional collocation methods for fractional differential equations (FDEs). It is also note¬ 
worthy that the choice of certain special JGL points with parameters related to the order of the 
equations can ease the implementation. We highlight that the use of the Bateman’s fractional 
integral formulas and fast transforms between Jacobi polynomials with different parameters, are 
essential for our algorithm development. 


1. Introduction 

Fractional differential equations have been found more realistic in modelling a variety of physical 
phenomena, engineering processes, biological systems and financial products, such as anomalous 
diffusion and non-exponential relaxation patterns, viscoelastic materials and among others. Typ¬ 
ically, such scenarios involve long-range temporal cumulative memory effects and/or long-range 
spatial interactions that can be more accurately described by fractional-order models (see, e.g., 

I3HG23031 USES and the references therein). 
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One challenge in numerical solutions of FDEs resides in that the underlying fractional integral 
and derivative operators are global in nature. Indeed, it is not surprising to see the finite dif¬ 
ference/finite element methods based on “local operations” leads to full and dense matrices (cf. 
[35 , l32l I3UI HU [15] .16, 42] [22] and the references therein), which are expensive to compute and 
invert. It is therefore of importance to construct fast solvers by carefully analysing the structures of 
the matrices (see, e.g., [431 131] 4. This should be in marked contrast with the situations when they 
are applied to differential equations of integer order derivatives. In this aspect, the spectral method 
using global basis functions appears to be well-suited for non-local problems. However, only limited 
efforts have been devoted to this very promising approach (see, e.g., [291 [301 EH EH m i), when 
compared with a large volume of literature on finite difference and finite element methods. 

Another more distinctive challenge in solving FDEs lies in that the intrinsic singular kernels 
of the fractional integral and derivative operators induce singular solutions and/or data. Just to 
mention a simple FDE involving RL fractional derivatives order p £ (0,1): (_?D/ t u)(x) = 1 for 
x £ (—1,1), such that u(— 1) = 0, whose solution behaves like u(x) ~ (1 + x) li . Accordingly, it only 
has a limited regularity in a usual Sobolev space, so the naive polynomial approximation has a poor 
convergence rate. Zayernouri and Karniadakis [49] proposed to approximate such singular solutions 
by Jacobi poly-fractonomials (JPFs), which were derived from eigenfunctions of a fractional Sturm- 
Liouville operator. Chen, Shen and Wang [9] modified the generalised Jacobi functions (GJFs) 
introduced earlier in Guo, Shen and Wang m, and rigorously derived the approximation results in 
weighted Sobolev spaces involving fractional derivatives. The JPFs turned out to be special cases 
of GJFs, and the GJF Petrov-spectral-Galerkin methods could achieve truly spectral convergence 
for some prototypical FDEs. We also refer to [45] for interesting attempts to characterise the 
regularity of solutions to some special FDEs by Besov spaces. It is also noteworthy that the 
analysis of spectral-Galerkin approximation in [29l I30j was under the function spaces and notion 
in |16j . and in [22], the finite-element method was analyzed for the case with smooth source term 
but singular solution. 

It is known that by pre-computing the pseudospectral differentiation matrices (PSDMs), the 
collocation method enjoys a “plug-and-play” function with simply replacing derivatives by PSDMs, 
so it has remarkable advantages in dealing with variable coefficients and nonlinear PDEs . However, 
the practicers are usually plagued with the dense, ill-conditioned linear systems, when compared 
with properly designed spectral-Galerkin approaches (see, e.g., ESS!). The “local” finite-element 
preconditioners (see, e.g., [25]) and “global” integration preconditioners (see, e.g., nil na [2oi in 
i46l [47] ) were developed to overcome the ill-conditioning of the linear systems. When it comes 
to FDEs, it is advantageous to use collocation methods, as the Galerkin approaches usually lead 
to full dense matrices as well. Recently, the development of collocation methods for FDEs has 
attracted much attention (see, e.g., [28] [50], :43, 17]). It was numerically testified in [28] [50] that 
for both Lagrange polynomial-based and JPF-based collocation methods, the condition number of 
the Caputo F-PSDM of order fj, behaves like 0(N 2/i ) which is consistent with the integer-order 
case. However, it seems very difficult to construct preconditioners from finite difference and finite 
elements as they own involve full and dense matrices and suffer from ill-conditioning. 

The main purpose of this paper is to construct integration preconditioners and new basis func¬ 
tions for well-conditioned fractional collocation methods from some suitably defined fractional 
Birkhoff polynomial interpolation problems. In [35], optimal integration preconditioners were de¬ 
vised for PSDMs of integer order, which allows for stable implementation of collocation schemes 
even for thousands of collocation points. Following the spirit of [46], we introduce suitable frac¬ 
tional Birkhoff interpolation problems at general JGL points with respect to both Caputo and 
(modified) Riemann-Liouville fractional derivatives (note: the RL fractional derivative is modified 
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by removing the singular factor so that it is well defined at every collocation point). As we will see, 
the extension is nontrivial and much more involved than the integer-order derivative case. Here, 
we restrict our attention to the polynomial approximation, though the ideas and techniques can be 
extended to JPF- and GJF-type basis functions. On the other hand, using a suitable mapping, we 
can transform the FDE (e.g., the aforementioned example) and approximate the smooth solution 
of the transformed equation, which is alternative to the direct use of JPF or GJF approximation 
to achieve spectral accuracy for certain special FDEs. 

We highlight the main contributions of this paper in order. 

• From the fractional Birkhoff interpolation, we derive new interpolation basis polynomials 
with remarkable properties: 

(i) It provides a stable way to compute the exact inverse of Caputo and (modified) 
Riemann-Liouville fractional PSDMs associated with “interior” JGL points. This 
offers integral preconditioners for fractional collocation schemes using Lagrange inter¬ 
polation basis polynomials. 

(ii) Using the new basis, the matrix of the highest fractional derivative in a collocation 
scheme is identity, and the F-PSDMs are not involved. More importantly, the resulted 
linear systems can be solved by an iterative method converging within a few iterations 
even for a very large number of collocation points. 

• We propose a compact and systematic way to compute Caputo and (modified) Riemann- 
Liouville F-PSDMs of any order at JGL points. In fact, we can show that the computation 
of F-PSDM of order k + p, with k £ N and /i £ (0,1) boils down to evaluating (i) F- 
PSDM of order p in the Caputo case, and (ii) a modified fractional integral matrix of 
order /i in the Riemann-Liouville case. Using the Bateman’s fractional integral formulas 
and the connection problem, i.e., the transform between Jacobi polynomials with different 
parameters, we obtain the explicit formulas of these matrices. 

The rest of the paper is organised as follows. The next section is for some preparations. In Sec¬ 
tion [3j we present algorithms for computing Caputo and (modified) Riemann-Louville F-PSDMs. 
In Sections |4][5l we introduce fractional Birkhoff polynomial interpolation and compute new basis 
functions. Then we are able to stably compute the inverse of F-PSDMs at “interior” JGL points 
and construct well-conditioned collocation schemes. The final section is for numerical results and 
concluding remarks. 


2. Preliminaries 

In this section, we make necessary preparations for subsequent discussions. More precisely, we 
first recall the definitions of fractional integrals and derivatives. We then collect some important 
properties of Jacobi polynomials and the related Jacobi-Gauss-Lobatto interpolation. We also 
highlight in this section the transform between Jacobi polynomials with different parameters, 
which is related to the so-called connection problem. 

2.1. Fractional integrals and derivatives. Let N and R be the sets of positive integers and 
real numbers, respectively, and denote by 

N 0 := {0} U N, R+ := {a <E R : a > 0}, R+:={0}UR + . (2.1) 

The definitions of fractional integrals and fractional derivatives in the Caputo and Riemann- 
Liouville sense can be found from many resources (see, e.g., [3S)[T2]): Nor p £ R + , the left-sided 
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and right-sided fractional integrals of order p are defined by 

(„«»)(*) = r^j l <•<“ >(I > = F(5) L i^h dv ' (2 ' 2) 

/or x £ (a, &), respectively, where T(-) zs i/ie Gamma function. 

Denote the ordinary derivative by D k = d k /dx k (with k £ N). In general, the fractional 
integral and ordinary derivative operators are not commutable, leading to two types of fractional 
derivatives: For p £ (k — 1, k) with fcgN, the left-sided Caputo fractional derivative of order p is 
defined by 

rx uW(y) 


u) {x) = a I k ~» ( D k u) (x) = j 


dy, 


(x - y)v~k+l 

and the left-sided Riemann-Liouville fractional derivative of order p defined by 


(*D!fu){x)=D k ( a I k -»u){x) = 


1 


j k r x 


u(y) 


-dy. 


(2.3) 

(2.4) 


T(k — p)dx k J a (x — yY k+1 
Note that if p = k £ N, we have „ ; D. k = = D k . 

Remark 2.1. Similarly, one can define the right-sided Caputo and Riemann-Liouville derivatives: 

(<?££«)(a) = (-1 ) k x I k b ~» (. D k u)(x ), (§D?u){x) = (-l) fc D k ( x I^u)(x). (2.5) 

With a change of variables: 

x = b + a — t, u(x) = v(a + b — x), x,t £ (a,b), 

one finds 

(tJ£v)(t) = ( a IPu){x), x,t£ (a, b), (2.6) 

and likewise for the fractional derivatives. In view of this, we restrict our discussions to the left¬ 
sided fractional integrals and derivatives. □ 


Recall that for p £ (k — 1, k) with k £ N, 


k ^ ( j) f \ 

(to? u)(x)= (aD£ u) (it) + Yi pTTT 




[x — a 


\3-^ 


(2.7) 


(see, e.g., mm, which implies 

(aD£u)(x) = (aD x u)(x), if u ij) (a) = 0, j = 0, - • - , fe — 1. (2.8) 

Moreover, there holds (see, e.g., [121 Thm. 2.14]): 

aDff a Ijf u(x) = u(x) a.e. in (a, b), p£ 1 + . (2.9) 

In addition, we have the explicit formulas (see, e.g., [I2J P. 49]): for real 77 > —1 and p £ R + , 


(x - a) v = + - aY+Y 


T(ri + p + l) 


and for p £ {k — 1 , k) with k £ N, 

f°’ 

aD(x -aY = l p( 7? + !) 


if T 7 e { 0 , ,k- 1 }, 


-(a;-a ) 17 72 , if 77 > k - 1 , 77 e K. 

,T{r]- p + 1) 

Similarly, for p £ (k — 1, k) with k £ N, and real 77 > —1, we have (cf. [38], P. 72]) 


( 2 . 10 ) 


( 2 . 11 ) 


(x - aY = 


rfa + i) 

T{r]- p+l) 


(x - aY~ k - 


( 2 . 12 ) 
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Hereafter, we restrict our attention to the interval A := (—1,1), and simply denote 

It := _i/£, c Dt := _$D£, R Dt = , x e A. (2.13) 

Apparently, the formulas and results can be extended to the general interval (a, b) straightforwardly. 

2.2. Jacobi polynomials and Jacobi-Gauss-Lobatto interpolation. Throughout this paper, 
the notation and normalization of Jacobi polynomials are in accordance with Szego pflj . 

For a, (3 € R, the Jacobi polynomials are defined by the hypergeometric function (cf. Szego [ill 
(4.21.2)]): 

Pi a ’ 0) {x) = “ITT— -^ 2 F 1 (-n,n + a + f3 + l;a + l- 1 ^—-x G A, n e N, (2.14) 
n!l(a + l) V 2 / 

and Pq C ‘’ i3 \x) = 1. Note that P^^\x ) is always a polynomial in x for all a, (3 gR, but not always 
of degree n. A reduction of the degree of Pn a ’^\x) occurs if and only if 


m:=-(n + a + (?)eN and 1 < m < n 
(cf. m P- 64] and .7]). Note that for a,/3 £ R, there hold 


P^\x) = P (^) (1) = r ^ r + Q + +| ) . 


(2.15) 


(2.16) 


For a, 13 > —1, the classical Jacobi polynomials are orthogonal with respect to the Jacobi weight 
function: o/“ ,/3 )(x) = (1 — x)“(l + x)^, namely, 


P^\x)P^\x)J^\x) dx = 7 < a ’«<W, 


(2.17) 


where 5 nn i is the Dirac Delta symbol, and 


M) = 2^+ 1 r(n + q + l)r(n + /3 + l) 

ln (2n + a + /3 + l)n! T(n + a + /3 + 1)' 1 ' 1 

However, the orthogonality does not carry over to the general case with a or 0 < —1 (see, e.g., 

[271 and [26] Ch. 3]). 

The following formulas derived from Bateman fractional integral formulas of Jacobi polynomials 
[4] (also see [21 P. 313], [41] P. 96] and [9]) are dispensable for the algorithm development. 

Theorem 2.1. Let p,s G R + , n £ No and x G A. Then for a € M and /3 > — 1, we have 

I p _{(l+xfP^\x)} = F( " + P + (1 + xf +p P^ +p \x), (2.19) 

T{n + p + p+l) 


and 


r D s _{{1 + x f+sp^-s,0+s) (g.)} = r( ^ + f + S+l \ l+xrP^(x). (2.20) 

1 (n + p + 1) 


As direct consequences of Theorem l2.il we have the following important special cases. 
Corollary 2.1. For a£l,p£ K + , n G No and x £ A, 

:(1 + x) p Pj L a ~ p ’ p \x); 


R D p _{(l+x) p Pi a ’ p \x)} = r(n + C +1 ) pi a+p ’0)(x). 


( 2 . 21 ) 


(2.22) 
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In particular, for p £ R + , n £ No and x £ A, 
I P -{Pu{x)} = 


n\ 


r(„ + „ + i 

R Di{(l+xYPp"\x)} = r( ’* + ^ 1| p,(4 


(2.23) 

(2.24) 


Remark 2.2. Remarkably, the formulas (I2.23I) - (I2.24I) link up the Legendre polynomials with the 
non-polynomials (l+x) p Pn~ p,p \x). They are referred to as the generalised Jacobi functions |19ll9.. 
and as the Jacobi poly-fractonomials [49] when 0 < p < 1. □ 

For a,/3> —1, let { Xj := x^j\cjj := be the set of Jacobi-Gauss-Lobatto (JGL) 

quadrature nodes and weights, where the nodes are zeros of (1 — x 2 )DP^’ l3 \x). Hereafter, we 
assume that {xj} are arranged in ascending order so that Xq = — 1 and Xn = 1. Moreover, to 
alleviate the burden of heavy notation, we sometimes drop the parameters a, f) in the notation, 
whenever it is clear from the context. 

The JGL quadrature enjoys the exactness (see, e.g., [39} Ch. 3]): 

,i n 

/ cj){x)uj^ a ^\x) dx = ^ 4>(xj)(jjj, V</> G P 2 N- 1 , (2.25) 

J ~ 1 3 =0 

where Vn is the set of all polynomials of degree at most N. Let X^u be the JGL Lagrange 
polynomial interplant of u £ C(A) defined by 


N 

(l N u)(x) = ^ u(xj)hj(x) £ Vn, 

3 =0 

where the interpolating basis polynomials {hj}jL 0 can be expressed by 

N 


(2.26) 


n —0 


with 


(«,« =7 («,/3) ( 0 < n < N - 1 




~(a,/3) 

Tn 


f — p(aj3)( \ 

L nj ~(ol,P) n 

In 

(2.27) 

2 +“ + £ +1 )7P>. 

(2.28) 


2.3. Transform between Jacobi polynomials with different parameters. Our efficient com¬ 
putation of fractional differentiation matrices and their inverses, relies on the transform between 
Jacobi expansions with different parameters. It is evident that for a, /3, a, b > — 1, 

V N = span {p&’fl : 0 <n<N}= span {p/“’ b) : 0 < l < N}. 

Given the Jacohi expansion coefficients of u £ Vn, find the coefficients {w| a,b ' ) } such that 

N N 

u{x) = £ u^P^Xx) = £ u\ a ' b) Pi a ' b) {x). (2.29) 

n —0 1=0 

This defines a connection problem (cf. [3]) resolved by the transform: 

= (a,P) C (a,b) (2.30) 

where and u l ' a ’ b ' > are column -(N +1) vectors of the coefficients, and j s the connec¬ 
tion matrix of the transform from {pX*’^} to One finds from the orthogonality (12.1711 

and (12.2911 that the entries of (. a ’h)c^ a ’ b \ i.e., the connection coefficients, are given by 


In ' (a,h) 

n 


p/ a ’ b) {x)P^ ] (s) W (a ’ b) (x)dx. 


(2.31) 
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Some remarks are in order. 

• By the orthogonality (12.1711 . we have = 0 for n < l, so the connection matrix is 

an upper triangular matrix. Therefore, (12.B0I) yields 


N 


u\ a ' b) = Y ( “ l/3) C, ( “’ b) 0 < l < N. 


(2.32) 


1=1 


In fact, we have the explicit formula of the connection coefficient (cf. [21 P. 357]) 

(a,0) c (a,t) ={2l + a + b+ 1) r(n + Q + l) r(Z + a + b+l) 
ln K T(n + a + /3 + 1) P(Z + a + 1) 


n—l 

E 

m—0 


(—l) m r(n + l + m + a + /3 + l)r(m + l + a + 1) 
m\(n — l — m)W(l + m + a + l)P(m + 21 + a + b + 2) 


(2.33) 


This exact formula is less useful in computation, as even in the Chebyshev-to-Legendre 
case, significant effort has to be made to analyze their behaviors and take care of the 
cancellations, when N is large (cf. [HE])- One can actually compute the connection 
coefficients by using the Jacobi-Gauss quadrature with (N + 1) nodes and with respect to 
the weight function o/ a > fe ). 

• In general, it requires 0(N 2 ) operations to carry out the matrix-vector product in (12.301) . In 
practice, several techniques have been proposed to speed up the transforms (see, e.g., [II [37, 
mm and the monograph [23] and the references therein). In particular, through exploiting 
the remarkable property that the columns of the connection matrix are eigenvectors of a 
certain structured quasi-separable matrix, fast and stable algorithms can be developed 
(cf. [23l [5j and the references therein). The interesting work [21] fully used the low-rank 
property of the connection matrix, and proposed fast algorithms based on rank structured 
matrix approximation. 


3. Fractional pseudospectral differentiation 

In this section, we extend the pseudospectral differentiation (PSD) process of integer order 
derivatives to the fractional context, and present efficient algorithms for computing the fractional 
pseudospectral differentiation matrix (F-PSDM). We show that 

(i) in the Caputo case, it suffices to evaluate Caputo F-PSDM of order /r G (0,1) to compute 
F-PSDM of any order (see Theorem 13.Ill : 

(ii) in the Riemann-Liouville case, it is necessary to modify the fractional derivative operator 
in order to absorb the singular fractional factor (see (13.81) 1. and the computation of the 
modified F-PSDM of any order boils down to computing a modified fractional integral 
matrix of order ^ G (0,1) (see Theorem 13.31) . 


3.1. Fractional pseudospectral differentiation process. It is known that the pseudospectral 
differentiation process is the heart of a collocation/pseudospectral method for PDEs (see, e.g., 
® ESI). Typically, for any u G Vn , the differentiation D k u is carried out via (12.261) in an exact 
manner, that is, 

N 

D k u(x) = Y u(xj)D k hj(x), k G N. 

7=0 


(3.1) 
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It is straightforward to extend this to the fractional pseudospectral differentiation. More precisely, 
for any u £ Vn, 

N 

(D^u){x) ='^2u{x j )D^h j (x), D* := p £ R+. (3.2) 

j=o 

However, in distinct contrast to (13.11) . we have D^u, D^hj ^ Vn, if p fL N. To provide some insights 
into this, we introduce the space: 

jff := {(l+s)> : V<t>eV N }, u£M, (3.3) 

and show the following properties. 

Lemma 3.1. For p £ (k — 1, k ) with k £ N, and for any u £ Vn, we have 

c Dfu £ F$lj? , R Dtu £ , (3.4) 

and 

c Dfu -4 0, R Dfu -A oo as i-i —1. (3.5) 

Proof. It is clear that 

D k u £ VN-k = span{P„ : 0 < n < N — k}. 

Thus, we derive from the definition (12.31) and (12.231) with p = k — p that 

c Df_u = I l f~ >J '(D k u) = (1 + ir) fc-/i <(>, for some <j>£VN-k- (3-6) 

Similarly, in the Riemann-Liouville case, we deduce from (12.231) that £ F^ ^. Then by the 

definition (12.41) . we obtain from a direct calculation that 

R Dfu = D k (I k ~^u) = (1 + cc) _M i/’; for some if £ Vn- (3-7) 

Thus, (13.41) is verified, from which (13.51) follows immediately. □ 

Remark 3.1. This implication of Lemma I5TT1 is that 

(i) if a FDE has a smooth solution, the source term might have a singular behaviour; 

(ii) conversely, for a FDE with smooth inputs, the solution might possess singularity. 

To achieve spectrally accurate approximation for some prototype FDEs pertaining to the latter 
case, the recent works [491 El proposed to approximate the singular solutions by using Jacobi 
polyfractonomials and general Jacobi functions, i.e., the basis of F^. □ 

Observe from (1331) that the Riemann-Liouville fractional derivative of any polynomial tends 
to infinity as xo —> — 1. This brings about some inconvenience for the computation of the related 
F-PSDM and implementation of the collocation scheme. This inspires us to multiply both sides 
of (13.21) by the singular factor (1 + x) M , leading to the modified Riemann-Liouville fractional 
pseudospectral differentiation: 

N 

( R Dfu){x) = ^u(x j )( R D , fh j )(x) where R Df := (1 + xy R D». (3.8) 

3=0 

With such a modification, we can recover the Riemann-Liouville fractional derivative values at 
Xi ^ -1 by 

( R Dfu){xi) = {l + x i )-> i ( R Dfu){x i ), 1 <i<N. (3.9) 

Correspondingly, we can define the modified factional integral and state some important properties 
as follows. 
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Lemma 3.2. Let u £ Vn and {hj} be the Lagrange interpolating basis polynomials at JGL points 
as before, and define 

It = (l + x)-' i Ilt, R Df := {1 + X y R Df, V/te R+. (3.10) 

Then we have 

Ifu. R D'fu £ Vn, oVm = span{/0/ij : 1 <j<N}, (3.11) 

where o Vn = {fi £ Vn '■ <(>(—1) = 0}. 

Proof. It is clear that by (13.711 . R Df_u £Vn, and by (12.2311 and (13.1011 . 

if{P n (x)} = —p(-^)(x), p £ M+. (3.12) 

Note that for any real p > 0, (x) is a polynomial of degree n (cf. [HI P. 64]). Thus, we 

have 

Vn = span {l^P n : 0 < n < IV} = span {ifhj : 0 < j < A7}, (3.13) 

and Ifu £ Vn- 

We now show (Ifhj)(— 1) = 0 for 1 < j < N. It is clear that {(1 + x)P}fi’ 1 ^ forms a basis 
of qVn, and by (12.191) with a = p and ft = 1, 

t{(l + x)P^ 1 \x)} = j^±^(l + x)P^\x). (3.14) 

Since hj £ o Vn, the identity (13.141) implies 1) = 0 for 1 < j < N. □ 

3.2. Caputo fractional pseudospectral differentiation matrices. As before, we use boldface 
uppercase (resp. lowercase) letters to denote matrices (resp. vectors), and simply denote the entries 
of a matrix A by Aij. Introduce the Caputo F-PSDM of order p : 

c D m £ rW+!)x(JV+ 1) ) C[)W = c Dfhj{xi). (3.15) 

In particular, for p = k S N, we denote D (k ^ = c D t ' k ' > and D = D 1 ' 1 ' 1 . 

Remarkably, the higher order Caputo fractional PSDM at JGL points can be computed by using 
the following recursive relation. 

Theorem 3.1. Let p £ (0,1). Then we have 

c D (k+n ) = c D {ti) D (k) = c D {n) D k ) k e N) ( 3 . 16 ) 

where D k stands for the product of k copies of the first-order PSDM at JGL points. 

Proof. For any u £ Vn, we have 

N 

u\x) = u {xi)h'i(x). (3.17) 

/—o 

Taking u = h' in (13.171) . leads to 

N 

tij( x ) =^2 h 'j(xi) h 'i{x), (3.18) 

i= o 

which, together with the definition (12.311 . implies 

N N 

C D^hj{x) = /i" M h"(x) = Y ti^x^iy^tifix) = Y hj(xi) c Dfhi(x). (3.19) 

1=0 1=0 

Taking x = Xi in the above, we obtain the matrix identity: 

c d G+u) = c D W D ^ ^ g ( 0i 


(3.20) 
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This leads to (13.161) with k = 1. Taking u = hj k \x) in (13.171) . we can derive the first identity in 
(13.161) in the same fashion. 

Using the property (see [39] Thm. 3.10]): 

D {k) =D k , fceN, (3.21) 

we obtain the second identity in (13.161) . □ 


It is seen from Theorem l3.1l that the computation of Caputo F-PSDM of any order at JGL points 
boils down to computing the first-order usual PSDM D (whose explicit formula can be found in 
e.g., [23]), and the Caputo F-PSDM °D^ with fx G (0,1). We present the formulas below. 


Theorem 3.2. Let {xj = }y^- 0 W ^ 1 a iP > ~ 1 and Xq = —1 be the JGL points, and 

let {wj = u}^’P} N _ q be the corresponding quadrature weights. Then the entries of C D < '^ with 
/i G (0,1) can be computed by 


c d[? = (1 + a*) 1 -" £ + Sli pfr 1 ' 1 -* (a *). 


(3.22) 


i=i 


for 0 < i, j < N, where 
1 W 


Slj 


= 2 51 (n + a + P + l) (a+1,/3+1) G ( -i 0) n-i tnj, tnj ■= ~^) p n a,f3) ), (3-23) 


1=1 — 1 


|(a+i,p+i) c -(o,o) n _j are Jacobi-to-Legendre connection coefficients, and are defined 

in (12.281) . In particular, for a = /3 = 0, we can alternatively compute the coefficients {s/j} by 


Slj = 


~{SjN + (-l)^jo -UjPt_ 1 (x j )}, 7z_i = 


7i-i 


21 - 1 


(3.24) 


To avoid the distraction from the main results, we provide the derivation of the formulas in 
Appendix [X] 


Remark 3.2. We see that in the Legendre case, we can bypass the connection problem. It 
is noteworthy that in [28], the Caputo F-PSDM of order n > 0 was computed largely by the 
derivative formula of P n and some recurrence relation of If_P„ built upon three-term recurrence 
formula of Legendre polynomials. As shown above, the use of the compact, explicit formula (12.231) 
leads to much concise representation and stable computation. □ 


3.3. Modified Riemann-Liouville fractional pseudospectral differentiation matrices. In¬ 
troduce the matrices: 

Rfilu) jlu) € R (JV+1)X(JV+1) where RfiM = ^d^hj)(xi), i[f = (Ifthj)(xi). (3.25) 
We can show the following important property similar to Theorem 13.II 


Theorem 3.3. Let {hj} be the JGL interpolating basis polynomials. Then for /x G (k — 1, k) with 
k £ N, we have 




where the entries of ^ are given by 
dk) 


D\j = (1 + *)"£>*{(! +x) k -»hj(x)}\ x=x ,, 0 < i,j < N. 


(3.26) 

(3.27) 
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Proof. By f|3. 1 1[) . we can write that for any u £ Vn , 

N 

(/y M zt)(x) = ^2(i k r ll u)(xi)hi(x), (3.28) 

1=0 

Multiplying both sides by (1 + a;) fe_/i , and using (13.101) . we find 

N 

{I k r*u){x) = 1 + x) k ~^hi{x), (3.29) 

1=0 

which implies 


N 

( R D^u)(x) =D k {I k r* l u)(x)=J2(i k T ,i u)(xi)D k {{l + x) k -»hi{x)}, (3.30) 

1=0 

for x £ (—1,1]. To remove the singularity at x = —1, we multiply both sides of (I3.30P by (1 + x ) M , 
and reformulate the resulted identity by the modified operator in (13.81) . leading to 

N 

( R D^ ) u)( X ) = J2(i-~ >1 U')(x l ){(l+xyD k {(l + x) k - ll h l (x)}y (3.31) 

1=0 

Taking u = hj and x = Xi in the above equation yields (13.261) . □ 


_ jjj) 

Observe from (13.26[) that it suffices to compute the modified fractional integral matrix I 
w (fc) 

with k — p £ (0,1), since D can be expressed in terms of the PSDM of integer order, e.g., for 
k = 1, 


-D U) = (1 - h)In+ i + A D, A = diag((l + x 0 ), • • • , (1 + x N )), 
where In+i is an identity matrix. 

Theorem 3.4. Let {xj = 0 with a,/3 > — 1 and Xq = —1 be the JGL points, and let 

{uij = be the corresponding quadrature weights. Then the entries of P ^ with p £ (0,1) 

can be computed by 


(3.32) 


W = E 


N 


l\ 


1 =0 


r (/ + /i +1) 


kjPi ^(.Xi), 0 <i,j<N, 


where 


N 


S, . _ V (.^)n(°<0)f . f . — p (a,P) (T A 

~ °2ra l ' n J, in 0 ■— ~(a,p) n b 


(3.33) 


(3.34) 


n=l 


7 n 


with {} b e i n g the Jacobi-to-Legendre connection coefficients, and defined in ()2.28D . 

In particular, if a = /3 = 0, we have sij = tij. 


Proof. It is essential to use the explicit formulas in Corollary 12.11 Accordingly, we expand the 
JGL Lagrange interpolating basis polynomials {hj} in terms of Legendre polynomials, and resort 
to the connection problem to transform between the bases as before. Equating (12.271) and the new 
expansion leads to 

N N 

hj (x) = J2 tnjP^'P) (x) = J2 0) - 0 < 3 < N \ 

n —0 1—0 


(3.35) 
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which defines a connection problem. Thus by (|2.32|) . 


N 


N 


% = E = E VM) 

n —0 n=l 

where we used the property: = 0 if n < l. Then it follows from (13.121) immediately that 

for n G (0,1), 

N 


4- =(^i)te) = £ 


v. 


Sij 


i=0 


r(/ + n +1) 


p[~^ (Xi), o <i,j<N. 


(3.37) 


This leads to the desired formulas. 

In the Legendre case, it is clear that the expansions in (13.351) are identical, so we have Sij = 


tij. 


□ 


We conclude this section by providing some numerical study of (discrete) eigenvalues of F- 
PSDMs. Observe from (13.51) that the first row of C D^ is entirely zero, so C D ^ is always singular. 
We therefore remove the “boundary” row/column, and define 




D := 


( c D ( m) ) 

V 1 ] n <i,j<N’ 


if M e (0,1), 


where C d[^ = ( c D l ^hj)(xi), (3.38) 

, .. , if II c II 91 

*3 )\<i,j<N-V 

which is invertible and allows for incorporating boundary condition(s). Similarly, we define 

>1 <i,j<N’ 


r D^ ] := 


if Me (1,2), 


(R Dij )i<u<N> if Me (0,1), _ (M) 


(xb^) 

V u )\ 


where R D i j = (^Dthj) (a:,). (3.39) 

’ij Ji<i,j<N-v if Me (1,2), 

In Figure 13.11 we illustrate the smallest and largest eigenvalues (in modulus) of these matrices. 
Observe that in both cases, the largest eigenvalue grows like 0{N 2fJ, ) : while the smallest one remains 
a constant in the Caputo case, and mildly decays with respect to N in the modified RL case. 




N N 

Figure 3.1. Maximum and minimum (in modulus) eigenvalues of F-PSDM with 
fi = 1.5. Left: Caputo. Right: modified Riemann-Liouville. 
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4. CAPUTO FRACTIONAL BlRKHOFF INTERPOLATION AND INVERSE F-PSDM 


As already mentioned, the condition number of the collocation system of a FDE of order p 
grows like 0(N 2/i ), so its solution suffers from severe round-off errors, and it also becomes rather 
prohibitive to solve the linear system by an iterative method. Following the spirit of [101146] . we in¬ 
troduce the Caputo fractional Birkhoff interpolation that generates a new interpolating polynomial 
basis with remarkable properties: 

(i) It provides a stable way to invert the Caputo F-PSDM in (13.3811 . leading to optimal frac¬ 
tional integration preconditioners for the ill-conditioned collocation schemes. 

(ii) It offers a basis for constructing well-conditioned collocation schemes. 

4.1. Caputo fractional Birkhoff interpolation. Let {a Oj = x^p}P 0 (with Xq = —1 and 

Xn = 1) be the JGL points as before. Consider the following two interpolating problems: 

(i) For p £ (0,1), the Caputo fractional Birkhoff interpolation is to find p £ Vn such that 

c Dfp(xj) = c Dfu(xj), 1 < j < N; p(- 1) = u(- 1), (4.1) 

for any u £ C[— 1,1] satisfying £ C(— 1,1]. 

(ii) For p £ (1, 2), the Caputo fractional Birkhoff interpolation is to find p £ Vn such that 

c Dfp(x j ) = c Dtu(x j ), l < j < N — 1] p(±l) = u(±l), (4.2) 

for any u £ C[— 1,1] satisfying < lD0u € C(—1,1). 

Remark 4.1. The usual Birkhoff interpolation is comprehensively studied in e.g., the monograph 
[33] . Typically, a polynomial Birkhoff interpolation requires at least one point at which the function 
and the derivative values are not interpolated consecutively. For example, consider a three-point 
interpolation problem: find p £ V 2 such that 

p(-l) = u(-l), p{0) = u{0), p(l) = u(l). 

It defines a Birkhoff interpolation problem, since the function value at x = 0 is not interpolated, 
as opposite to the Hermite interpolation. Due to the involvement of Caputo fractional derivatives, 
we call m and (14.211 the Caputo fractional Birkhoff interpolation problems. □ 


As with the Lagrange interpolation, we search for a nodal basis to represent the interpolating 
polynomial p. More precisely, we look for Qj £ Vn such that 
(i) for p £ (0,1), 

C Df Qj{xi) = Sij, l<i<iV; Qj{—1) = 0, 1 <J<N, (4.3) 

with Qq = 1; 


(ii) for p £ (1,2), 


c DfQ^x i ) = 6 lj , 1 < i < N — 1; Q?(± 1) = 0, 1 < j < N - 1, (4.4) 

with Qq{x) = (1 — x)/2 and Qf^ix) = (1 + x)/2. 

Then, we can express the Caputo fractional Birkhoff interpolating polynomial p of SU and s 
respectively, as 

N 

P(x) =u{-l) + '^2 c Dfu(x j )Q‘*(x), p £ (0,1), (4.5) 

i =1 


and 


N -1 


V{x) = 1) + ‘b-uix^Qfix) + ^p-u{ 1), p £ (1,2). 


3 =1 


2 


2 


(4.6) 
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Therefore, {Qj} are dubbed as the Caputo fractional Birkhoff interpolating basis polynomials of 
order p. 

Introduce the matrices 


Q (M) = 


(Qif) 


ij >i<l,j<N’ 

)Wi 


if p e (0,1), 


where Q[f =Qf(xi). 


if MS (12) 

Remarkably, the matrix Q ^ is the inverse of c d[^ defined in (13.3811 . 


(4.7) 


Theorem 4.1. For p £ (k — 1, k) with k = 1, 2, we have 

q(F c d= c d (F) Qf/x) = l N+1 _ k , 


(4.8) 


where I^+i-k is the identity matrix of order N + 1 — k. 


Proof. We just prove (|4.8D with p £ (0,1), as the case p £ (1, 2) can be shown in a similar fashion. 
Since Qj £ Vn and Qj{— 1) = 0 for 1 < j < N, we can write 

N N 

Qj( x ) — Qj( x i)hi(x) = 5 ~2Qj(xi)hi(x ), 1 <j<N, 

1=0 1=1 

where {hi} are the Lagrange interpolating basis polynomials associated with JGL points. Thus, 

N 

c Df Q^x) = Qj( x i)° D - H x )- 

1=1 

Taking x = X{ for 1 < i < N in the above equation, we obtain (14.811 with p £ (0, 1) from (14.311 
straightforwardly. □ 


4.2. Computing the new basis {Qj}- The following property plays a crucial role in computing 
the new basis {Qj}, which follows from Lemma f3. II 

Lemma 4.1. Let be the JGL points with = —1 and xi v = 1. Then for p £ (k — 1, k ) 

with k = 1,2, we have 

D k Q^x) = R D k S»{Q^) k ~*hj{x)), l<j<N+\-k, (4.9) 

where {hj}^l([ 1— are the Lagrange-Gauss interpolating basis polynomials associated with the JGL 
points {xj}™Z[ ~ k , that is, 

hj £ Vn-u , hji x i) = Sij for 1 < i, j < N + 1 — k. (4-10) 

Proof. Since Qj £ Vn, we obtain from Lemma [34] that c Df Q 1 } £ ■ Noting that {fi J }jZ 1 1 ^ k 

forms a basis of Vn-r, so by m , 

N+l-k 

c DfQ^x)= J2 cij {l+x) k -^hi{x), l<j<N+l — k. 

1=1 

Letting x = x% and using the interpolating conditions, we find that cij = (1 + xi)^' Thus, we 
obtain 

°Dt<? j {x)= (±±lLy-\{x), l<j<N+l — k. (4.11) 

By the definition (12.311 . we have c D l f = I k V^ 1 D k , so using (EH), we obtain (14.91) from (14.111) 
immediately. □ 
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With the aid of (14.91) . we are able to derive the explicit formulas for computing the new basis. 
We provide the derivation in Appendix iBl 

Theorem 4.2. Let {xj = w ( a >P > ~ 1 and xq = —Xn = —1) be the 

JGL quadrature nodes and weights. Then {Qj} can be computed by 
(i) For p£ (0,1), 


N-l 


QIW = 


_ 1 y 

C 1 + Xj ) 1 ^ ^0 


r(i-n + 2) 
i\ 


itj / Pi{x)dx, 1 <j<N, 


where 


N-l 


4 = E 


, i f Cj - pTE-i) 

7^1 /3 + 1 P^\x 3 ) 




with Cj = 1 for 1 < j < N — 1, and cjv = a + 1. 
(ii) Por G (1,2), 


w = 


1 -g r( '-; +3) 4i., M , i<i<jv-i, 


(1 + ^ 

v ;=o 


where 


N—2 


4 = E ( ^E ( r 2 ’ 2_M Ep 


T=l 


i r (^ - i)P^ i/3) (-i) 


^ 7 ^ ) \2(^ + l)P^fe) 


P^-«(-lVo 


and 


{l + x 3 )P { «’ P) (l) 
2(a + l)P^\ Xj ) 


Ji a ^(l)w N + Pi a ^( Xj )u 


3)^3 ft 


^i(. x ) J (t-l)Pi(t)dt + j (x-t)Pi(t) dt. 

Here, { ( “.«C' I ( £ _fc,fc_, * ) } are f/ie connection coefficients as defined in (12.4111 . 


(4.12) 

(4.13) 

(4.14) 

(4.15) 

(4.16) 

(4.17) 

(4.18) 


Remark 4.2. Observe from (14.1311 and (14.1611 that if we take (a, /3) = (p — k, k — p) with fc = 1,2, 
we have = £y, so (14.121) and (14.1511 have the simplest form. Thus, it is preferable to choose 
these special parameters. □ 


5. MODIFED RL FRACTIONAL BlRKHOFF INTERPOLATION AND INVERSE F-PSDM 

We introduce in this section the fractional Birkhoff interpolation involving modified Riemann- 
Liouville (RL) fractional derivatives which offers new polynomial bases for well-conditioned col¬ 
location methods for solving FDEs with Riemann-Liouville fractional derivatives. Moreover, we 
are able to stably compute the inverse matrix of r d\^ defined in (13.3911 . However, this process 
appears more involved than the Caputo case in particular for p G (1, 2). 
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5.1. Modified Riemann-Liouville fractional Birkhoff interpolation. Like the Caputo case, 
we consider the modified Riemann-Liouville fractional Birkhoff interpolating problems (i)-(ii) as 
defined in m-m with r D in place of c Df. Similarly, we look for the interpolating basis 
polynomials {Qj}™ = 0 C Vn such that 
(i) for /t G (0,1), 

R DHQS(xi)= 0, 1 <i<N; Qq(—1) = 1, 


R DfQ^{x i ) = 6 ij , 1 <i<N, Qf(- 1) = 0, 
(ii) for n G (1,2), 

R D- Qo(xi) = 0, 1 < * < iV — 1; Qq(—1) = 1, 

R DtQ$(x i )=6 ij , l<i<N-l, Q;(± 1) = 0, 
R DtQ^( Xi ) = 0, 1 < i < N — 1] Q%{- 1) = 0 

Then for any u E Vn, we can write 

N 


l<j<N ; 

QS( 1) = 0, 

1 < j < N - 1, 

Q5v(i) = i- 


(5.1) 


(5.2) 


u(x) = u(—l)Qft(x) + y ^ R D^u(xj)Qt(x) (for/r G (0,1)) 


i =i 
IV-1 


(5.3) 


= tt(-l)Qft(a;) + Y R D R u{xj)Q R (x) +u(l)Q%(x) (for p G (1,2)). 
i=i 

Introduce the matrices generated from the new basis: 


~(u) 

Q = 


(«“) 
.(«“) 




if V G (0,1), 
if At G (1,2), 


;(/*) 


where =Q^{xi). 


(5.4) 


Like Theorem ITT! we can claim that ^ is the inverse of r d\^. As the proof of the theorem 
below is very similar to that of Theorem 14.11 we omit it. 

Theorem 5.1. For p G (k — 1, fc) with k = 1,2, we have 


Q 


(F R ^(F 


^in — Mn b? — JjV+l-fe) 


(5.5) 


where I^+i-k is the identity matrix of order N + 1 — k. 


5.2. Computing the new basis {Q^} ._ Q . The following lemma is very useful for the computa¬ 
tion, whose proof is provided in Appendix [Cl 


Lemma 5.1. Let p G (fc — 1, k) with k = 1,2. Then for any f G o'PiV) the fractional equation 

R Dfu(x) = f(x), u(- 1) = 0, (5.6) 

has a unique solution u G oVn of the form 

u(x) = If{{l+x)~^f(x)}. (5.7) 

In particular, for any u G Vn, we have 

R Dfu(— 1) = 0 if and only if u(— 1) = 0. (5-8) 

For clarity of presentation, we deal with two cases: (i) p G (0,1) and (ii) p G (1,2), separately. 
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5.2.1. {Qj}^- o with h G (0,1). Using the properties (13.111) and (15.81) . we obtain from the interpo¬ 
lating conditions in (15.11) that 

( R &LQfi(x) = h j (x), 1 <j<N; ( R D»_QZ)(x)=Zh 0 (x), (5.9) 

where {hj} are the JGL interpolating basis polynomials defined in ( 12 . 271 ) . and £ is a constant to 
be determined by Qq(— 1) = 1 . Note that thanks to ( 15 . 81 ) . the condition Qj(—1) = 0 is built-in, as 
hj{— 1) = 0 for 1 < j < N. We summarise below the explicit representation of the new basis. Once 
again, we put the proof in Appendix iDl 


Theorem 5.2. Let {xj = x^f\oJj = (with a, ft > —1 and xq = —1) be the JGL 

quadrature points and weights. Then with \x G (0,1) can be computed by 


Q?(s) = 0E r( * H + With iu = (5.10) 

1=0 n=l 

where Co = l/r(l-/i), & = lforl<j<N, {^C^} are the connection coefficients defined 
in Subsection 12.31 and 

Vj (5-11) 


N 


CUj 

tn i = J 


7« 


with 7 ^“’^ being defined in (12.281) . 


Remark 5.1. If (a, ft) = (p,—p) with p G (0,1), we have tij = tij, so has the simplest 
form. □ 


5.2.2. {Q1}% with /i G (1,2). It is essential to derive the identities like (15.91) . Indeed, using 
(13.111) and (15.81) . we obtain from the interpolating conditions in (15.21) that 

( R DfQZ)(x) = (T Q + K0 x)h 0 (x), Qo(-l) = 1, Qg(l) = 0; (5.12) 

( R D?_Q*){x) = ±±p-h j {x), Qj(l) = 0, 1<J<N-1; (5.13) 

( R DfQ%)(x) = T N (l+x)h 0 (x), Qjy(l) = 1, (5.14) 

where { hj are the Lagrange interpolating basis polynomials at JGL points {TCjjy^ 1 , that is, 
hj{x) G 'Pat-i, hfixi) = Sij, 0 < i, j < N — 1. (5.15) 


In (15.12l) - (|5.14l) . {rj}jL 0 and kq are constants to be determined by the corresponding conditions 
at x = ±1, e.g., Qj(l) = 0 in (15.131) . It is noteworthy that thanks to (15.81) . the interpolating 
condition: Qj(— 1) = 0 is built in ( R Df Qj)(— 1) = 0 for 1 < j < TV. 

In what follows, we shall use the three-term recurrence relation of Jacobi polynomials (cf. 55) 
(3.110)]): 


xP^-^ix) = a l+1 P^-^\x) + hP^-^ix) + ci-xP^-^ix), l > 0 , 

where c_i = 0, fi G (1, 2), and 

l -h 2 1 — 2 fl (/ + /i)(Z — fJj H - 1) 

ai+1 = 27T3’ bl = (21 + l)(2l + 3) ’ Ci_1 = (l + l)(2l + l) ■ 


(5.16) 

(5.17) 


As before, it is necessary to expand {hj} in terms of Jacobi polynomials with different parameters 
by using the notion of connection problems, so as to use compact and closed-form formulas to 
compute the new basis. We state below the connections of three expansions, and postpone the 
derivations in Appendix [El 
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Lemma 5.2. Let , ojj = u N,j } -o ( w ^ d a,/3 > —1 and xq = — Xn = —1) be the JGL 


1 3- 

quadrature nodes and weights, and let be the Lagrange interpolating basis polynomials 

c j} N ~o 1 ^ 1 r ' h 

JV-1 N- 1 


associated with {xj}P q 1 defined in (15.1511 . Then for p G (1,2), we have 


hj(x) = Y Qnj PY 0 \x) = Y Qlj p i ^ 0 ) 


n —0 


1= 0 


N—2 


(5.18) 


= soj + Y ( X + rfPi 11 ’ 1 M \ x )> 0 < j < AT - 1, 


z=o 


where goo = 1 and goj = 0 for 1 < j < N — 1. Moreover, the coefficients can be computed by 

1 




7n 

AT-1 


-1, 


eij = E Qnj , 0 < l, j < N - 1, 

n=l 


(5.19) 

(5.20) 

(5.21) 


and by the backward recurrence relation: 


a 1 ^ bj, ~ (~ 1 A Ci A 

Qij — Qij Z?i+i,j Qi+2,j, i = N — 3, N — 2, ■ ■ • , 1, 

0-2 0.2 

„ 1 „ 1 6jV_ 2 + 1 „ 

QN—\,j = - QN-l,j, QN-2,j = - QN—2,j - QN-l,j, 


(5.22) 


Oat-1 

where {ai,bi,Ci} are given in (15.171) . 


O-N-2 


OJV-2 


With the above preparations, we are ready to derive the explicit formulas of the new basis 
{Qj}j - 0 h G (1,2). We refer to Appedix [F] for the derivation. 

Theorem 5.3. Let { Qij,Qij } be the coefficients defined in Lemma 15.21 and denote 


d i := r ^ + + 2 i)!^ ’ : = (1 + ir)i^ (0 ’ 1) (ai). 


(5.23) 


Then at JGL points with fi G (1,2) can be computed by 

(i) for j = 0, 


N -1 


N—2 


Qo( x ) =1 + ( r 0 - E d iSioTi{x) + p^ 1 ^^ E ^fft+i.owW 


where tq — 


1 


z=o 

1 


N—2 


r (i -m) 12 r(i- M )^ 


, AT-1 


(5.24) 


E <a + i,o f / E 


z=o 


(ii) /or j = 77, 


AT-1 


AT-1 


QnOa) = ^ E d i dw / E d i ; 


;=o 


1=0 


(5.25) 
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(ii) for 1 < j < N - 1, 


N—2 


Qj( x ) = «+ij Vi( x ) + 


Xj + To 
3 3 1=0 


Xj + Tj 


N-2 


55 d isi+ij <1 l+ , 2 o ^ n-iw+iW + ^m{x) 


;=o 


Z + 2 


l + 1 

l + 1 - /i 


Q-iV^-iOc) ?, 


(5.26) 


where {ai,bi,ci} (with c_i = 0) are defined in (15.171) . and 


iN—2 

tj = -1 + /ur(2 - /x)f?ij / 55 d f 
' i- o 


(5.27) 


Remark 5.2. We see from (15.211) that if (a,/3) = (fi, 1 — //) with /r G (1, 2), the connections 
coefficients are not involved, so Qj has simpler form. □ 


6. Well-conditioned collocation schemes and numerical results 

In this section, we apply the tools developed in previous sections to construct well-conditioned 
collocation schemes for initial-valued or boundary-valued FDEs, and provide ample numerical 
results to show the accuracy and stability of the methods. 

6.1. Initial-valued Caputo FDEs. To fix the idea, we first consider the Caputo FDE of order 
MG (0,1): 

C D0n(x) + A (x)u(x) = f(x), x £ (— 1,1]; u(— 1) = u_, (6.1) 

where A, / are given continuous functions, and it_ is a given constant. The collocation scheme is 
to find ujv G TV such that 

c Dfu N (xj) + A (xj)u N (xj) = f(xj), l<j<N; u N (- 1) = u— (6.2) 

The corresponding linear system under the Lagrange basis polynomials {h :j } (L-COL) becomes 

( C D^+A)u = f, (6.3) 

where C D[^ is defined as in (13.381) . A = diag(A(xi), • • • , A(xjv)), and 

u= (u N (x i),--- ,u N (x N ))\ f = (f(x i) - U- C D l fh 0 (x 1 ), ■ ■ ■ ,f{x N ) - U- C D l fh 0 (x N )y. (6.4) 

The collocation system under the Birkhoff interpolation basis polynomials {Qj} in (14.31) (B- 
COL) becomes 

(l N + AQ^)v = g, (6.5) 

where 

N 

u N (x) = U- + 55' y l ( 5l ( x )’ v = (vir ■ ■ , v n) , (6.6) 

2=1 

and g = (/(xi) — it_A(xi),--- , /(xjv) — rt_A(xjv)) t - It is noteworthy that different from (16.31) . 
the unknowns of (16.51) are not the approximation of u at the collocation points, but of the Caputo 
fractional derivative values in view of (14.51) . 

Thanks to Theorem 14.II we can precondition (16.3|) and obtain the PL-COL system: 

(l N + Q^A)u = Q^f. 


(6.7) 
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In the computation, we take X(x) = 2 + sin(25x) and u(x) = i(—2(1 + xY) with /j, = 0.8 in 

m, where the Mittag-Leffler function is defined by 


(- 2 ') 


oo 


E 

71=0 


Y(na + /3)" 


( 6 . 8 ) 


In view of Remark PI we choose the JGL points with (a, ft) = (/it — 1,1 — n) = (—0.2,0.2). We 
compare the condition numbers, number of iterations (using BiCGSTAB in Matlab) and conver¬ 
gence behaviour (in discrete L 2 -norm on fine equally-spaced grids) of three schemes (see Figure 
16.11) . Observe from Figure 16.11 (left) that the condition number of usual L-COL divided by N 2fl 
behaves like a constant, while that of PL-COL and B-COL remains a constant even for N up to 
2000. As a result, the latter two schemes only require about 8 iterations to converge, while the 
usual L-COL scheme requires much more iterations with a degradation of accuracy as depicted in 
Figure Ed] (middle). We record the convergence history of three methods in Figure [6J] (right), and 
observe that two new schemes are stable even for very large N. 




Figure 6.1. Comparison of condition numbers (left), iteration numbers against 
errors (middle), and errors against N at convergence in log-log scale (right) for 

(ED- 


6.2. Boundary-valued Caputo FDEs. We now turn to the boundary value problem: 

c D^u(x) + X\(x) c D l !_u(x) + X 2 (x)u(x) = f(x), x £ (—1,1); 

(6.9) 

u{— 1) = U-, u(l) = u + , 0 < v < /j,, /i e (1,2), 

where Ai, A 2 and / are given continuous functions, and u± are given constants. 

With a pre-computation of the Caputo fractional differentiation matrices of order /r and v in 
Subsection we can formulate the L-COL scheme as (16.31) straightforwardly. The counterpart 
of (16.51) i.e., the B-COL scheme, can be formulated as follows: find (cf. (14.61) 1 

N ~i l-x 1 _)— x 

u N (x) = u* N (x) + ^2 Vj Qj{x), (1,2); u* N (x) :=—^u-+—^u+, (6.10) 

3 =1 


such that 


(Jv-i+AiQ H +A 2 Q (m) )« = 9 


(6.11) 
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Figure 6.2. Comparison of condition numbers (left), iteration numbers against 
errors (middle), and errors against N at convergence in log-log scale (right) for 

s 


where A* = diag^fyi), • • • , A^ar^-i)), * = 1,2, Q^ = c D u _Q^{xi), 1 < i,j < N - 1, and 
g = (f(x i) — g*(*i), • • • , f(xN-i) ~ q*(xN—i)y with g* = \i C D v _u* n + A 2 fyjfy Note that the entries 
of Q ^ can be evaluated by Theorem 12. 11 (12.71) and (I4.15I) - (I4.18I) . Here, we omit the details. 


Remark 6.1. If Ai = 0 and A 2 is a constant, we can follow (46J Proposition 3.5] to justify the 
coefficient matrix of (16.Ill) is well-conditioned. Indeed, thanks to Theorem 14.11 the eigenvalues a 
of I iv —1 — A 2 Q (/i) satisfy 


1 + A 2 A 


-1 

max 


< a < 1 + A 2 A 


-1 

min’ 


where A max and A m i n are respectively the largest and smallest eigenvalues of • Since A m i n = 
0(1) (see Figure EH] (right)), the condition number of Jjv-i — is independent of N. □ 


Like (16.71) . we can precondition the L-COL scheme by Q ^ which leads to the PL-COL system. 
In the following comparison, we set g = 1.9, ^ = 0.7 and (a, f3) = (—0.1, 0.1) (cf. Remark I4~121) . 
and take 

Ai(x) = 2 + sin(47ra;), A 2 (x) = 2 + cosx, (6.12) 

and 

u{x) = e 1+x + (1 + x) 6+4/7 - 2(1 + z) 5+4/7 , (6.13) 

where we can use the formula 


°D»_e 1+x = (1 + x) k ~^E hk+1 _ tl { 1 + x), g€ ( k - 1, k), k = 1,2, 
to work out f{x). 

Once again, we observe from Figure [6721 that the new schemes: B-COL and PL-COL are well- 
conditioned, attain the expected convergence order about 10 iterations, and lead to stable compu¬ 
tation for large N. 


6.3. Riemann-Liouville FDEs. Consider the Riemann-Liouville version of (16.91) : 

R Dtu(x) + \\(x) R D v _u(x) + \ 2 (x)u(x) = f(x), x e (—1,1); 
u(— 1) = U-, u(l) = u+, 0 < v < g, g S (1, 2), 

where Ai, A 2 and / are given continuous functions, and u± are given constants. 


(6.14) 
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For a better treatment of the singularity, we consider the modified Riemann-Liouville fractional 
collocation scheme: find un G Vn such that 

R D^u N (xj) + Xi(x j ) R D'iu N (xj) + A 2 (xj)u N (xj) = f(xj), 1 < j < N; 
u N (-l) = u-, ujv(1) = u+, 0 < v < n, /t G (1, 2), 


(6.15) 


where Ai = (1 + xY~ v \i, A 2 = (1 + x) M A 2 , and / = (1 + x)^f. 

Here, we just focus on the collocation system using the new basis in (15.31) . that is, 

N-l 

u N {x) = u_Q£(x) + u + Q%{x) + ^2 v j Qj(x), /x G ( 1,2). (6.16) 

3 =1 

Then one can write down the B-COL system in a fashion very similar to (16.111) with only a change of 

~ (v) ^ (fl) 

basis. Correspondingly, we denote the matrix of the linear system by A := In-i+AiQ +A 2 Q 
We first show that the B-COL scheme enjoys spectral accuracy (i.e., exponential convergence), 
when the underlying solution is sufficiently smooth. For this purpose, we take 

u(x) = e" (1+x) - -e“ 2 ^-^, (6.17) 

and Ai, A 2 to be the same as in (16.121) . In Figure [6731 we plot discrete L 2 -errors for various pairs 
of (/x, v ) of the B-COL schemes for both Caputo and Riemann-Liouville fractional boundary value 
problems (BVPs) (16.91) and (16.141) under the same setting. We observe the exponential decay (i.e., 
0(e~ cN ) for some c > 0) of the errors. Both schemes take about 10 iterations to converge, while 
much more iterations are needed and severe round-off errors are induced if one uses the standard 
L-COL approach. 




Figure 6.3. Errors against N of the B-COL schemes for Caputo fractional BVP 
(16.91) (left) and Riemann-Liouville fractional BVP (16.141) (right) with various /r, v, 
where Ai, A 2 are given in (16.121) . and the exact solution u is given in (16.171) . 


We further test the new B-COL method on (16.141) with smooth coefficients but large derivative: 

Ar(i) = 1 + e" 1000 * 2 , A 2 (*) = 1 + e _1000 ( x + 0 - 2 ) 2 , (6.18) 

and with the exact solution having finite regularity in the usual Sobolev space: 

u{x) = ( - (1 + xY/2) + 1 ~ E ^- 2,t ') (i + x )-l. 


2 


(6.19) 
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We tabulate in Table 16.11 the discrete L 2 -errors, number of iterations and the second smallest and 
largest eigenvalues (in modulus). Once again, the scheme converges within a few iterations even 

for very large N. In fact, as we observed from Figure EH] (right), the smallest eigenvalue of R D in 
in (13.391) still mildly depends on N. As a result, the condition number of A grows mildly with 
respect to N. However, it is interesting to find that the eigenvalues in modulus of A (denoted by 
which are arranged in ascending order) are concentrated in the sense that 

0(1) = M < \<Jj\ < \<?n-2\ = 0(1), 2 < j < iV — 2. (6.20) 

Thanks to this remarkable property, the iterative solver for the modified Riemann-Liouville system 
is actually as fast as the previous Caputo system where the coefficient matrix is well-conditioned. 

Table 6.1. Errors, number of iterations and concentration of eigenvalues of A 


N 


/i = 1.5, v = 

0.6 


= 1.9, v = 

0.7 

1^2 | 

|o\N 2 

Iters 

Errors 

Wi 1 

|o\N 2 

Iters 

Errors 

8 

0.7 

1.0 

7 

8.58e-03 

0.4 

1.0 

6 

2.66e-03 

16 

0.6 

1.1 

12 

2.03e-03 

0.2 

1.7 

8 

3.69e-04 

32 

0.6 

1.2 

12 

5.39e-04 

0.2 

3.0 

8 

5.49e-05 

64 

0.6 

1.3 

12 

1.31e-04 

0.1 

5.0 

8 

7.46e-06 

128 

0.6 

1.5 

12 

3.24e-05 

0.1 

7.4 

8 

1.06e-06 

256 

0.6 

1.6 

12 

8.08e-06 

0.1 

9.9 

8 

1.53e-07 

512 

0.6 

1.7 

13 

2.02e-06 

0.1 

12.7 

8 

2.21e-08 

1024 

0.6 

1.7 

13 

5.04e-07 

0.1 

21.1 

8 

3.69e-09 


6.4. Concluding remarks. In this paper, we provided an explicit and compact means for comput¬ 
ing Caputo and modified Riemann-Liouville F-PSDMs of any order, and introduced new fractional 
collocation schemes using fractional Birkhoff interpolation basis functions. We showed that the 
new approaches significantly outperformed the standard collocation approximation using Lagrange 
interpolation basis. 

As a final remark, we point out two topics worthy of future investigation along this line, which we 
wish to explore in forthcoming papers. The first is to analyze the fractional Birkhoff interpolation 
errors and understand the approximability of the new interpolation basis functions from theoretical 
perspective. The second is to extend the idea and techniques in this paper to study the fractional 
collocation methods using the nodal basis {(l+x) M /ij(x)/(l-|-a;j) /i } (see [50] . i.e., the counterpart 
of Jacobi poly-fractonomials [35] and generalised Jacobi functions [45]). 


Appendix A. Proof of Theorem 13.21 
We expand {h'A in terms of Legendre polynomials, and look for {szj} such that 


N 


N 


h , J (x) = -Y,(ri + a + p + l)t nj Pt + 1 1 ’ fi+1 \x) = J2^Pi-i(x), 0 < j < N, (A.l) 


1 = 1 


(A.2) 


where {t n j} are given in (12.271) and we used the derivative formula (cf. j411 ): 

DP^\x) = l( n + a + p + 1 )Pti 1 ' 0+1 \x), n > 1. 

By (I2.29I) - (I2.32I) . we find that 

N N 

s ij = ^E( B+a+ ^ 1 ) ( “ +W1 * C f’J-i t ”i = ) J2 {n+a+ /3+l) ( “ +1 ’^ +1) C'|°’ 1 °^_ 1 t n j , (A.3) 

n= 1 n=l — 1 
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for 1 < l < N and 0 < j < N. In view of (|A.1|) . we can use the first formula in (|2.23|) to derive 


N 




1 = 1 


N 


= (i + ^) 1 - M E 


!)■ D (M-i.i-A') 


(A.4) 


' S H p i 


i -1 


{Xi) . 


r(z +1 — p) 

This ends the derivation of (I3.22l) - (13.23l) . 

We now derive (13.241) for the LGL case. Using the orthogonality of Legendre polynomials, 
integration by parts and the exactness of LGL quadrature (cf. (12.251) '). we obtain from (IA.1I) that 

N 


1/1 if 1^ V 1 

sij = - / h'Ax)Pi-\{x) dx = - hj(x)Pi-i(x)\ - E /u^OA-iO^M 

7i-i 7-i 72-i l iZ o J 


72-1 


UjJV + ( 1) ^jPl-liXj)}, 


2-1 


(A.5) 


where we used the properties: hj(xi ) = Sij and P;_i(±l) = (±1) 

Appendix B. Proof of Theorem 14.21 


Since hj £ pN-k , we can write 


N-k 


N-k 


Kj(z) = E Zn3 P n a ’ P \ x ) = E hP^-^ (*), l<j<N + l-k. 

n =0 1=0 

As before, if one can work out {£ nj }, then by (|2.29[) - (|2.32|) . 

N-k 


£ _ \ ^ (cKj/SVAM - k,k — fl) £ 


(B.l) 


(B.2) 


n=l 


As to be shown later, inserting (IB. II) into (14.91) . we can derive from (12.241) with p = k — p the 
desired formulas. Thus, it remains to find {£„,,} in (IB. II) . We proceed separately for two cases. 

(i) For p £ (0,1), we obtain from the orthogonality (12.171) . the exactness of JGL quadrature 
(12.251 ) . and the interpolating condition (14.101) that 

,1 . N 


£nj ~ 


Aa,p) 

An 


^ h j {x)P^\x)uj { - a ^\x)dx = -^'Yshj{x i )p( a ' l)) {xi)u i 

An ’ i=0 (B.3) 


In 


{hj{-l)P^\- l)wo + P^HxjPj}, 1 < j < N. 


Now, we evaluate hj(— 1). Since {A,} are associated with the interpolating points {x.j}jL ,, which 
are zeros of (1 — x)DP^ L ’^\x), we have the representation: 

h 3 (x) = 


(1 -x)DP^ ) (x) 


1 <j<N. 


(x - Xj)D{{l - x)DP^’ 0) (x)}\ x _ " 

Recall the Sturm-Liouville equation of Jacobi polynomials (cf. [H] (4.2.1)]): 
-(1 -x 2 )D 2 P ( *' P \x) = {P-a-{a + P + 2)x}DP { *’ 0 \x) + \ 

where A^*’^ = N(N + a + /? + 1). It follows from (IB. 51) that 

A a P)[ 1\ _ \ (a./3) r>(a,P)/ ,1 \noW) 


(«>/9) p(“./3) 
N r N 


2(/3 + l)DFr Pj (-l) = -X ( n 0 ) Pn’ 0) 


-{l-x 2 )D 2 P^\ Xj )=\ 


A 7V r N 


(-1), 2(a + 1)TP;' W (1) = A 
\xj), 1 <j < N -1. 


(x), 

( 1 ), 


(B.4) 


(B.5) 


(B.6) 
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Using the property: (1 — Xj)DP^’^ (xj) = 0 and (IB. 61) . we compute from (1B.4|) that 


p(a,t 3)/_,\ 

hj{- 1 ) =--—^ i<j<N, 


(B.7) 


where c 3 = 1 for 1 < j < N — 1, and cjv = a + 1. Substituting (IB. 71) into (IB. 31) yields (14.141) . 
Inserting (IB. II) into (14.91) , we derive from (12.241) with p = 1 — p that 


N-l 


DQj(x) = 


1 -^' r(i -; + 2) {„P, W , l<i<JV. 


0- + X3) 1 11 


1=0 


(B.8) 


In view of Qj(— 1) = 0, a direct integration of (IB.8I) leads to (14.121) . 

(ii) For p £ (1, 2), (IB.31) reads 

Znj = -^{h 3 {-l)P^X-l)u 0 + ^(l)Pi“’«(lViv + l<j<N-l. (B.9) 

7n 

We need to evaluate hj{± 1). Note that in this case, {h 3 } are associated with the interior JGL 
points ; which are zeros of DPj^’^Xx), so we have 


Hj(x) = 


DP^\x) 


{x — Xj)D' 2 P^’^ ) (xj) 




(B.10) 


Thus using (IB.5I) - (IB.6I) leads to 

M-i) = - 


1 -Xj p^X- i) 


2(/5 + l) P^' P Xxj) 


, h 3 { 1) = - 


l + x 3 P^X 1) 

2(« + l) P^Xxj) 


(B.ll) 


Substituting (IB.Ill) into (IB. 91) yields (14.171) . 

Similar to case (i), inserting 4EH> into (|4.9 p , we derive from (|2.24|) with p = 2 — p that 


N—2 


D2 « w =(i +Ij) ^ 


l—y 

r ,\2—/i / 


Z=0 


T(l-p + 3) 

l\ 


iijPXx). 


(B.12) 


Solving this equation with the boundary conditions: Qfj 1 (±1) = 0, we obtain in (14.181) and the 
desired formula (14.151) . 


Appendix C. Proof of Lemma I5TT1 

We carry out the proof by directly verifying that u{x) in (ED is the desired polynomial solution. 
It is evident that for any / £ qPn, we can write 


N—1 

f{x)=Y,In{l+x)P^Xx), (C.l) 

n =0 

where the coefficients {/„} are uniquely determined. Using (12.191) with p = p, a = p and /3 = 1 — p, 
leads to 

u(x) =I^{(l+x)-»f(x)} = Y)f ^^( 1 + :r ) P » 0 ’ 1) ( a '’)’ 

n =0 ' ' 


(C.2) 
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which implies u £ qVn- Recall that R D^_ = (1 + x) pR D p . Thus, acting R D P on both sides of (1C.21) . 
we obtain from (12.2011 and (1C.Ill immediately that 


N -1 


R D»U(X) = E r( i ^ / ) fn R D^{(l+x)P^\ X )} 


(n + 1)! 

= Y J fn{l + x)P^ ) {x)=f{x). 


n —0 
N—l 


(C.3) 


n—0 


Therefore, u(x ) in (EZD verifies m ■ The uniqueness follows from R D^u(x) = 0 implying u(x) = 0. 

We now turn to (15.81) . The above verification shows that if /(—1) = 0, i.e., R D p _u{— 1) = 0, 
then u(— 1) = 0. Hence, it suffices to show if u(— 1) = 0, then R D p _u{— 1) = 0. For this purpose, we 
expand 

JV-l 

u{x ) = E Rn + X ) P n 0,1) ( X )’ (C.4) 

n —0 

where {u n } can be uniquely determined. Like the derivation of (1C.31) . acting R D P and using (12.2011 . 
we obtain 

~ N ~ X (n 4- 1 ’ll 

R D P _u{ x) = E p( n + 2-n) ^ + ( c - 5 ) 


n—0 


which implies R D p _u(—l) = 0. 


Appendix D. Proof of Theorem 15.21 
We intend to use the compact identity deduced from (12.201) . that is, 


R Dt{P n (*)} = 


T(n — n + 1) 


P^~ p \x), n> 0, (0,1). 


(D.l) 


This inspires us to expand {hj} (resp. Q p ) in terms of {P/ M ’ E (resp. {P/}). Following (13.3511 
(13.36[h we have 


N 


N 


h j( x ) = E ^ P ^’ = E ^’^In M) tnj , 


1=0 


a=l 


and 


N 


Qj (a) = E ( x )> 0<j< N. 


(D.2) 


(D.3) 


z=o 

Inserting (1D.2I) - (ID.Hl) into (15.911 . we obtain from (ID.Ill immediately that for 0 < l < N, 

tu — n + i) ~ „ . Ar ^ ru-n + i)^ 

Qij =-1 < 3 < N, qio =- — -£ t w . 

m = qs 

Using (15.81) . the formula (12. 121) . and definition (13.101) . we obtain from (15.91) that 


(D.l) 

Thus, it remains to determine the constant £. Setting Qq(x) = Qq(x) — 1, we have Qq(—1) = 0. 


o = ( R 3 p Qo)(—i) = £ - ( R nft i) | x=1 = e - 


i 


F(1-m)’ 


so £ = 


1 


r (l -M)‘ 


(D.5) 


This ends the proof. 



































FRACTIONAL COLLOCATION METHODS 


27 


Appendix E. Proof of Lemma 15.2 


We first derive the coefficients in (15.19l) - (15.20ll . By the orthogonality (12.1711 . the exactness of 
JGL quadrature (12.2511 . and the interpolating condition (15.151) . we have 


Qnj = -T^ay f h j {x)P^\x)u)^\x)dx = —Y h 3 { Xl )P^ p) (®j)wi 

'Y™ ’ J-l 7" ’ i=0 


(E.l) 


(«,/ 3 ) 

In 


{P^\x 3 )u 3 + h 3 (l)P^\ 1) WJV }, 0 < n, j < JV — 1. 


Since {fty} are associated with the JGL points {p, ; which are zeros of (1 + x)DP^’^\x), we 

have the representation: 


hj{x) = 


(1 + x)DP ( £' P \x) 


{x-x 3 )D{( 1 + x)DP^\x)}\ x=x 
A direct calculation using (IB. 61) leads to 


0<j<N-l. 


(E.2) 


^ 0 (1) = JJvUL, L(1) = —— Y ^ , 1 < j < iV — 1. (E.3) 

Thus, we obtain lj5.19M5.20l) by inserting them into (IE. 111 . 

Thanks to 

N-l IV-1 

hj(x) = Y Qnj = Y Sij P^’^^ix), o < j < N - 1, 


(E.4) 


n —0 


1=0 


we solve the connection problem and obtain from (12.32I) - (I2.31I) the formula (15.2111 . 

It remains to derive (15.221) . Applying the three-term recurrence relation (15.161) to the last 
expansion in (15.181) . we obtain the connection 


TQj — Qji Qj — {Qoji " ■ ■ > QN—i,j) i Qj — (Qoji j QN—i,j ) j 


(E.5) 


where T is an upper triangular matrix with only nonzero entries on diagonal and two upper 
diagonals: 

Too = 1, Tu = dj; Ti t i + i = bi + 1; 2 = Cj. (E.6) 

Solving the linear system by backward substitution leads to (15.221) . 


Appendix F. Proof of Theorem 15.31 

We first use Lemma IQl to solve (I5.12D - (I5.14D and find the expressions of the constants therein. 
It’s more convenient to reformulate (15.121) as: find Qq(x ) = Q q(x) + 1 such that 


(R D ^_ 


Qo)(x) = ( 


= T o- 


(1 + x)ho(x) + 


r(i-f*)' 

where we used (I2.12|) . (13.81) and (15.81) to derive 


P(l-M) 


R Dt 1 = 


1 


t(i-m) 

Using Lemma l5Tl and (12.101) . we obtain 
1 


, K 0 = To 


(ho(x) — 1), Q£( 1) = -1, (F.l) 


(F.2) 


r(i -nY 


Qo(x) = ( 


= To - 


r(i 


)^-{(l + x) 1 ^ho(x)} + p( 1 1 _^ / -{( 1 + a: ) ^Cho(x) - 1 )}. (F.3) 
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As Qo(l) = — 1) we have 


r(i ~ h ) t 0 = i - 


^{(1 + x)-^ho(x) - 1 )} | a=1 + r(i - n) 


i^{{i + xy-vh 0 (x)}\ x=1 

Following the same argument, we derive 

Qn( x ) = t n I^{(1 + x) 1 ~ IJ -h 0 {x )} 1 t n = 

and for 1 < j < N — 1, 


7^{(l + x) 1 -/*/jo(x)}| a . =1 ’ 


Qj(x) = ^ 2^. (rjl-Ul + x) ^hj(x)}+!-{( 1 + x) ^xhj 


Xj + Tj 


T-i = - 


I-{( l + x) ^xhj{x)}\ x=1 
11{(l+x)-»hj(x)}\ x=1 


(F.4) 

(F.5) 

(F.6) 

(F.7) 


We now evaluate fractional integrals of hj. Using the last two expansions with j = 0 in (15.18[) . 
and the identity (12.191) with p = p,a = p and /? = 1 — p, we obtain 

^-{(1 +x) 1 ~ fl ho(x)} = Y r ^ / + + 2 ^| /i ' > &o (! + x)P l i ° A \x), 


;=o 


7^{(l + x) »(h 0 (x) - 1)} = Y r< ~// + | 2 i \\^ g»+i.o (! + x ) p i 0,1) ( x )- 


(F.8) 


;=o 


(* + !)! 


Noting that Pi 0,1 ^(l) = 1 (cf. [41]), we obtain from (IF. 41) and (IF. 81) the value of r 0 in (15.241) . and 
the expression of Q q (x) follows from (IF. 31) immediately. 

Similarly, we obtain from (IF. 51) and (IF. 81) the expression of Q^{x) in (15.251) . 

We now turn to Qj{x) with 1 < j < N — 1. Once again, using (12.191) (with p = p, a = p and 
/3 = 1 — p) and (15.181) . leads to 

^-{(1 +x)~^hj(x)} = (1 + x) Y r ^ + + 2 1 )| M ' > &+i,j p i 0,1) ( x )> 


1=0 
N-2 . 


A{(1 + *T'AM}U = 2 •£ 


(F.8) 


where we used p/ 0,1 ^(l) = 1. Moreover, by (12.191) . (15.181) and (15.161) - (15.171) . 

N-2 

I-{0 + x)~ ll xhj{x)} = Yj Ql+i,j I-{0 + a;) 1_A ‘xP/' i ’ 1_ ^(x)} = (1 + x)x 


^ T(l + 2- p) 
— 77 i iu — 


1=0 


(1 + 1)! 


1=0 

j l -b 2 — fi 

l 1 + 2 


y+ b,pp»+J-ti-c-tpv? w>, 


(F.10) 


l + 1 - /i 

where c_i = 0. Using the property: p[^\l) = 1 and (15.171) . we find from a direct calculation and 
(IF.91) that 

^{(l+xy^xhjix)}]^ = 2{l-p)T(2-p)eij+2Y r ^ + + 2 1 )i^ ft+ij n . 

= 70 {(1 + x)hj (x)} | x=1 - 2p r (2 - p) g 1:i . 

Inserting (IF.9I) - (IF.11D into (IF.6I) - (IF.7D . we derive the formulas (I5.26I) - (I5.27I) . 
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